
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <omp.h>
//////////////////////////////
#include <cuda.h>
#include <cuda_runtime.h>
//////////////////////////////

#define REPEAT 256

// Function that runs on GPU
__device__ __inline__ float f(float x,float y){
  return  sinf(cosf(sinf(sqrtf(powf(cosf(x*x+y*y),1.5f)))));
}


// Kernel: apply the function to all elements of an array
__global__ void apply_function_REPEAT(float * __restrict__ a,
				      const float * __restrict__ b){
  int id   = blockIdx.x * blockDim.x + threadIdx.x;
  for (int count=0;count<REPEAT;count++) a[id] = f(a[id],b[id]);
}


int main(void){
  long unsigned int N = 1024L*1024*256;
  float *a = (float *)malloc(N*sizeof(float));
  float *b = (float *)malloc(N*sizeof(float));
  struct timespec ts0,ts1,tsi;
  double time;
  
  fprintf(stderr,"*** GPU ***\n\nFilling arrays\n");
  for (long unsigned int i=0;i<N;i++){
    a[i] = (float)rand()/(float)(RAND_MAX);
    b[i] = (float)rand()/(float)(RAND_MAX);
  }

  clock_gettime(CLOCK_REALTIME,&ts0);
  fprintf(stderr,"Processing arrays\n");

  
  // Allocate arrays on GPU
  float *d_a1,*d_b1,*d_a2,*d_b2;
  cudaSetDevice(0);
  cudaMalloc(&d_a1,(N/2)*sizeof(float));
  cudaMalloc(&d_b1,(N/2)*sizeof(float));
  cudaSetDevice(1);
  cudaMalloc(&d_a2,(N/2)*sizeof(float));
  cudaMalloc(&d_b2,(N/2)*sizeof(float));
  clock_gettime(CLOCK_REALTIME,&tsi);
  time = (tsi.tv_sec-ts0.tv_sec + (tsi.tv_nsec-ts0.tv_nsec)/1e9);
  fprintf(stderr,"\tAllocation       : %.2f\n",time);
  
  // Transfer initial data to GPUs
  cudaSetDevice(0);
  cudaMemcpy(d_a1,a,(N/2)*sizeof(float),cudaMemcpyHostToDevice);
  cudaMemcpy(d_b1,b,(N/2)*sizeof(float),cudaMemcpyHostToDevice);
  cudaSetDevice(1);
  cudaMemcpy(d_a2,a+N/2,(N/2)*sizeof(float),cudaMemcpyHostToDevice);
  cudaMemcpy(d_b2,b+N/2,(N/2)*sizeof(float),cudaMemcpyHostToDevice);
  clock_gettime(CLOCK_REALTIME,&ts1);
  time = (ts1.tv_sec-tsi.tv_sec + (ts1.tv_nsec-tsi.tv_nsec)/1e9);
  fprintf(stderr,"\tCopy HOST -> GPU : %.2f\n",time);

  // Perform processing on GPU
  cudaSetDevice(0);
  apply_function_REPEAT<<<N/2/128,128>>>(d_a1,d_b1);
  cudaSetDevice(1);
  apply_function_REPEAT<<<N/2/128,128>>>(d_a2,d_b2);
  cudaSetDevice(0);
  cudaDeviceSynchronize();
  cudaSetDevice(1);
  cudaDeviceSynchronize();
  clock_gettime(CLOCK_REALTIME,&tsi);
  time = (tsi.tv_sec-ts1.tv_sec + (tsi.tv_nsec-ts1.tv_nsec)/1e9);
  fprintf(stderr,"\tProcessing       : %.2f\n",time);
 
  // Retrieve the result from GPU
  cudaSetDevice(0);
  cudaMemcpy(a,d_a1,(N/2)*sizeof(float),cudaMemcpyDeviceToHost);
  cudaSetDevice(1);
  cudaMemcpy(a+N/2,d_a2,(N/2)*sizeof(float),cudaMemcpyDeviceToHost);
  cudaSetDevice(0);
  cudaDeviceSynchronize();
  cudaSetDevice(1);
  cudaDeviceSynchronize();
  clock_gettime(CLOCK_REALTIME,&ts1);
  time = (ts1.tv_sec-tsi.tv_sec + (ts1.tv_nsec-tsi.tv_nsec)/1e9);
  fprintf(stderr,"\tCopy GPU -> HOST : %.2f\n",time);

  time = (ts1.tv_sec-ts0.tv_sec + (ts1.tv_nsec-ts0.tv_nsec)/1e9);
  fprintf(stderr,"Total : %.2f seconds\n\n\n",time);

  exit(0);
}
